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ABSTRACT 

We estimate the conditions for detectability of two planets in a 2/1 mean-motion resonance 
from radial velocity data, as a function of their masses, number of observations and the signal- 
to-noise ratio. Even for a data set of the order of 100 observations and standard deviations of 
the order of a few meters per second, we find that Jovian-size resonant planets are difficult to 
detect if the masses of the planets differ by a factor larger than ~ 4. This is consistent with the 
present population of real exosystems in the 2/1 commensurability, most of which have resonant 
pairs with similar minimum masses, and could indicate that many other resonant systems exist, 
but are presently beyond the detectability limit. 

Furthermore, we analyze the error distribution in masses and orbital elements of orbital 
fits from synthetic data sets for resonant planets in the 2/1 commensurability. For various mass 
ratios and number of data points we find that the eccentricity of the outer planet is systematically 
over estimated, although the inner planet's eccentricity suffers a much smaller effect. If the initial 
conditions correspond to small amplitude oscillations around stable apsidal corotation resonances 
(ACR), the amplitudes estimated from the orbital fits are biased toward larger amplitudes, in 
accordance to results found in real resonant extrasolar systems. 

Subject headings: celestial mechanics, planets and satellites: general 



1. Introduction 

Although the population of extrasolar plan- 
ets continues to increase rapidly, the number of 
multiple-planet systems with members in mean- 
motion resonances (MMRs) shows a much slower 
growth rate. The existence of a 3/1 commen- 
surability in 55Cnc has been recently questioned 
by the five planet fit of Fischer et al. (2008) 
whose best solution now indicates a non-resonant 
motion for 55Cnc-c and 55Cnc-d. Even in the 
2/1 MMR, the most dynamically important and 
populated commensurability, the number of con- 
firmed planetary systems is still restricted to four 
well known cases: GJ876, HD82943, HD73526 and 
HD128311. Gozdziewski et al. (2007) proposed 
that HD160691 may also have two planets in the 



2/1 MMR, although this is unconfirmed, and just 
two years before the best fit solution seemed to 
favor a 5/1 MMR between the same two plan- 
ets (Gozdziewski et al. 2005). Fig. [1] shows the 
distribution of exoplanet pairs according to their 
mass ratio and orbital period ratio. Except for 
the vicinity of the 2/1, they appear more or less 
at random. 

However, this picture may change in the near 
future. Recently, two possible new resonant sys- 
tems have been proposed based on radial velocity 
RV observations from the Geneva group. Laskar 
and Correia (2009) have found a planetary sys- 
tem around HD60532 which appears to be trapped 
in a 3/1 MMR, while Correia et al. (2009) show 
a similar behavior for the orbital solutions of a 
pair of planets around HD45364, this time in 
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around a 3/2 commcnsurability. Given the volatil- 
ity of orbital fits for near-resonant configurations 
(e.g. 55Cnc) it is perhaps too early to treat these 
cases as confirmed. Nevertheless, their impor- 
tance is unquestionable, specially since they would 
populate commensurabilities which are presently 
empty. 

For several years a large proportion of resonant 
planets was expected as a consequence of an as- 
sumed large-scale planetary migration of exoplan- 
ets due to interactions with the gaseous disk. In 
fact, the existence of resonant systems has many 
times been advanced as observational evidence 
that such a migration actually took place, and 
many planets were formed farther from the star 
than their present location. Hydrodynamical sim- 
ulations seem to indicate that resonance trapping 
in low-order commensurabilities (particularly the 
2/1 MMR) arc high- probability outcomes for a 
wide spectrum of planetary masses, initial condi- 
tions and disk parameters. Several works, particu- 
larly focused on GJ876 (see, e.g., Kley et al. 2005), 
point that the present configuration of planets b 
and c can only be explained assuming such a sce- 
nario. 

Three explanations have been presented re- 
cently to account for the lack of a larger resonant 
population. One possibility is that not all plan- 
ets approaching the 2/1 resonance may have been 
captured. If the mass ratio between the outer and 
inner planet was sufficiently small (of the order of 
the Saturn over Jupiter ratio), then a very fast or- 
bital decay of the outer smaller body (e.g., Type 
III migration) may have avoided resonance cap- 
ture in the 2/1 and led to a later trapping in the 
3/2 (Massct & Snellgrove 2001). A similar effect 
has been proposed by Morbidclli & Crida (2007) as 
a first step to explain the current orbital architec- 
ture of the outer planets of our own solar system. 
However, it is not clear under what circumstances 
such a fast migration would occur (see D'Angelo 
& Lubow 2008). Alternatively, turbulence effects 
in the gaseous disk may have caused significant 
orbital perturbations in the decaying bodies to in- 
hibit resonance trapping (Adams et al. 2008). 

A second possibility may lie in the survival rate 
of resonant planets during their evolution within 
the commensurability. As shown originally by Lee 
and Pealc (2002), once inside the resonance do- 
main, tidal interactions with a gaseous disk will 



increase the eccentricities of the planetary bodies 
until a disruptive collision or ejection of one body 
occurs. This is due to the topology of the stable 
ACR families in the eccentricity domain (see, e.g., 
Beauge et al. 2006, Hadjidemetriou 2008) within 
the 2/1 MMR. The only way two resonant planets 
may survive a large scale orbital migration is if the 
driving mechanism introduced a significant damp- 
ing of the orbital eccentricities, leading to equilib- 
rium values of these elements comparable with the 
observed values (Kley et al. 2005, Beauge et al. 
2006) . Although this effect is expected from nom- 
inal disc parameters, especially if an inner disk is 
assumed (Crida et al. 2008), there is no evidence 
that this must be true in all cases, and perhaps 
most of the systems originally trapped in the 2/1 
could have been ejected. Moorhead and Adams 
(2005) proposed such a mechanism to explain the 
present semimajor axis and eccentricity distribu- 
tion of exoplanets. 

A third possibility, little considered up to now, 
is that the apparent lack of resonant planets may 
be due to detectability limitations. Recall that 
resonant motion causes an almost periodic repe- 
tition of the RV curve of the star which, under 
certain circumstances, may not allow a good sep- 
aration of components. Perhaps, many additional 
systems may actually lie within the 2/1 resonance 
but are currently not discernible due to limited 
observations or small signal-to-noise ratios. Re- 
cently, Anglada-Escude et al. (2008) discussed 
cases where a two-planet resonant system with al- 
most circular orbits may appear masked as a sin- 
gle planet in an elliptical orbit of eccentricity e. 
However, this effect appears to be possible only 
in cases where the outer resonant planet is much 
more massive than its inner companion and the 
single-planet solution has a low eccentricity (e.g., 
e < 0.2). In this paper we address this question 
in more detail constructing a detectability crite- 
rion valid for any mass ratio and not restricted to 
quasi-circular orbits. 

We also discuss the errors in the planetary 
masses and orbital parameters of those resonant 
systems which can actually be detected. Although 
this is different problem, it shares many common 
points and can be studied using the same ap- 
proach. Of the four systems currently inhabit- 
ing the 2/1 MMR, only GJ876 has a dynamically 
stable best fit, while the others are characterized 
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Table 1: Orbital Parameter for Three ACR Solu- 
tions in the 2/1 MMR for Different Mass Ratios 
mijmx. Masses are in units of Jovian mass and 
angles in degrees. The central mass is mo = IMq. 



Fig. 1. — Distribution of pairs of consecutive plan- 
ets in multiple-planetary systems, according to 
mass ratio over ratio of orbital periods. MMRs are 
indicated by the vertical red lines, whose length is 
inversely proportional to the order of the MMR. 
The data correspond to most current orbital fits. 

by orbits leading to a disruption of the system in 
time scales much smaller than the age of the star. 
Although stable orbital solutions are possible for 
these troublesome cases, having rms similar to the 
best fit, they usually correspond to large ampli- 
tude ACR, not easily compatible with a smooth 
orbital migration (Sandor et al. 2007, Crida et al. 
2008). 

2. Orbital Fits with Synthetic Data Sets 

In order to avoid the problems of estimating er- 
rors in observational data with unknown solution, 
we will work with synthetic data sets of RVs. We 
will assume the existence of two planets of masses 
mi and mi orbiting a star mo with semimajor axes 
<Zi, eccentricities e^, mean longitudes \ and lon- 
gitudes of pericenter vdi. The index i — 1 will be 
used for the inner planet, while i = 2 will denote 
the outer body (i.e. a\ < 02). We assume that 
both bodies share the same orbital plane oriented 
edge-on with respect to the observer. 

From the nominal solution, we generate a syn- 
thetic RV curve mimicking the star movement 
around the barycenter of the system. This curve 
is the sum of two periodic signals, each with am- 
plitude Kj related to the j'-th planet. Our RV 
data set will only cover a few orbital periods, and 



we will assume no significant effects from mutual 
gravitational perturbations. 

Once the synthetic curve is generated, we con- 
struct a discrete sampling choosing N observation 
times ti distributed randomly according to a ho- 
mogeneous distribution. At each point, we calcu- 
late a RV value as a random displacement of the 
nominal V r (ti); this displacement follows a Gaus- 
sian distribution with constant variance a 2 . The 
resulting synthetic data set will consist of three 
columns {ti,V ri ,a) and will be used as input in 
our orbital fitting procedure (Beauge et al. 2008). 

2.1. Two-Planet Systems in a 2/1 Mean- 
Motion Resonance 

Consider two planets in a small-amplitude oscil- 
lation around a stable ACR in the 2/1 MMR. The 
question is how the uncertainties in the orbital fit 
will affect the observed motion of the system. As 
examples, we have chosen three different configu- 
rations, each corresponding to a stable ACR with 
different mass ratios m2/mi. The first two corre- 
spond to ACR of type (0, 0), while the last displays 
an asymmetric corotational behavior. Masses and 
initial orbital elements are given in Table [TJ In all 
cases the initial conditions lead to a small ampli- 
tude oscillation (approximately 5 degrees) around 
the respective ACR. 

Typical examples of synthetic RV data sets with 
cr = 10 m s^ 1 are shown in Fig. [2] Although all the 
generated data sets are constructed with the same 
standard deviation cr, the signal-to-noise ratio ap- 
pears different for each mass ratio. To understand 
this effect, let us write the total amplitude of the 
RV curve as amp(K-) = K\ + K-2- Expressing Ki 
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Fig. 2. — Black continuous line shows synthetic 
RV curves for the three ACR presented in Table 
[TJ Mass ratio decreases from top to bottom. Open 
circles are N — 165 fictitious points with random 
time distribution and standard deviation a = 10 
m s . 

in terms of the mass and orbital elements (see, 
e.g., Beauge et al. 2007) for an edge-on coplanar 
system, we can write (up to second order in the 
masses, and for quasi-circular orbits): 

amp (lO = Ki + K 2 = —n\a\ + — n 2 a 2 , (1) 
m m 

where rij are the mean motions. For planets in 
the vicinity of a 2/1 MMR, we can simplify this 
expression as: 

an*(V P )*(l+y 5 — J— max, (2) 

where {n 2 /iii){a 2 / a\) ~ y/l/2. The quantity in- 
side the brackets represents the increase in the 



RV amplitude due to the presence of the outer 
(resonant) planet. This term tends to unity for 
m 2 — > 0, and shows a linear dependence with the 
mass ratio. Thus, larger values oim 2 jm\ will pro- 
duce a larger RV signal and, assuming a fixed ob- 
servational standard deviation a, will result in a 
larger signal-to-noise ratio. 

Equation ^ is also a rough indication of which 
planet dominates the RV signal. The critical mass 
ratio is given by m 2 /m\ ~ \/2 — 1.26. For smaller 
values, the RV amplitude of the inner planet mi is 
larger, and the RV curve appears as a perturbed 
signal with primary period equal to the orbital pe- 
riod of the inner planet (i.e., 2-k /n{). One example 
is shown in the bottom frame of Fig. [21 Conversely, 
if the mass ratio is larger than ~ 1.26, the signal 
from the outer planet becomes larger and the dom- 
inant period in the RV curve is given by 2n/n 2 (see 
top graph of Fig. [J). The middle frame represents 
a transition region in which both components are 
of similar magnitude. 

3. Detectability Criteria 

Cumming (2004) presented a simple procedure 
to estimate the detectability of single-planetary 
systems, given the number of data points N and 
K/a ratio, as a function of the desired false-alarm 
probability (FAP). Although the original formu- 
lation was developed for circular orbits and for 
large data sets, it serves as a starting point for the 
extensions shown below. Throughout this paper 
we will assume that the observational time frame 
covers at least one orbital period of the planetary 
masses. 

Imagine that we have two fits of a given RV 
data set (V*) with N points. Each fit is assumed 
to contain v = N — M degrees of freedom (i.e., 
M free parameters) and to yield a squared sum of 
residuals 

N 

Q = J2(V i -W(t i )) 2 (3) 
l 

where W(t) is the adopted modefj]. As usual, the 
Vi are assumed to be statistically independent nor- 
mal variates V* = N(W(U),a 2 ). The variance a 



We prefer to use Q instead of \ 2 ( = Q/ a2 ) because of 
the ambiguity introduced on this notation by its improper 
usage in astronomy where x 2 is often " normalized" and not 
the same quantity used in statistics. 
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(unknown) is time- independent but, if necessary, 
the definition of Q can be changed to introduce 
weights and take into account different variances 
(see Ferraz-Mello, 1981). 

We adopt the subscripts a and b to the first and 
second fit, respectively. In addition, we assume 
that the model a is embedded into model b (M a < 
Mb). The improvement of the goodness of fit may 
be characterized by the power z defined as (see 
Cumming 2004) 

Vb Qa — Qb 



Qb 



(4) 



This formula extends to the general case the 
" floating- mean power" introduced in the study of 
periodograms by Cumming et al. (1999). When 
dealing with a Gaussian white noise, the probabil- 
ity distribution functions of this statistic is given 
by the Fisher-Snedecor distribution F Ua - VbtVb (z). 

The FAP (or probability of error) is the prob- 
ability of getting one result just by chance when 
working with a white noise. When the result is 
the solution of the minimum problem Q — Q m in, 
we must keep in mind that it does not refer to a 
random choice of the parameters, but to the max- 
imum z obtained from a given number of trials. 
Following Cumming (2004), the FAP T(z) of the 
best fit result z is given by 



T{z)=M J F Va ^ b!Ub (z)dz, 



(5) 



where M. is roughly the number of independent 
trials. In the case of a periodogram constructed 
with equidistant data points, the periodogram is 
a Fourier transform and A4 = N/2. For more 
involved problems, there is no simple expression 
and this quantity is often determined from Monte 
Carlo simulations. When the data are unequally 
spaced, empirical formulae are used. We men- 
tion the rules introduced by Cumming (2004) and 
by Quast (Ferraz-Mello and Quast, 1987). In 
the cases studied in this paper, we have found 
that M. ~ N seems to give a better agreement 
with Monte Carlo simulations. This is the value 
adopted through the present work. Equation ([5]) 
can then be inverted to give the detectability limit 
Zd corresponding to a user specified FAP^, i.e., 
solving 

FAP d = T{z d ) (6) 
to obtain the necessary Zd- 



3.1. Synthetic random samples 

In order to study how the detectability varies 
with the parameters and the quality of the avail- 
able measurements, we may study a simple model 
b in which we assume that the measurements are 
given by Vi — Wb(ti) + e; where Wb(t) is a chosen 
model (circular, Keplerian with e ^ 0, etc.) and 
is a Gaussian noise (normal variate) iV(0, of) un- 
corrected with Wb- As model a, we use the mean 
of the data, assumed to have been averaged to 
zero beforehand. We may use orthonormal func- 
tions (see the Appendix) to get an approximation 
for the statistic z corresponding to the best-fit so- 
lution^. We can rewrite equation ([50)1 as 

t = ^(\m\+(N-irt_ 1 ) i (7) 

Q b J 
which may be simplified to read 

Z! = ^ 8 

Q b 

where we have introduced the auxiliary statistic 



Z\ = z + 



v a - v h 



and the constant 

Vb 



A = 



N(W? 



+ N- 1 



(9) 



(10) 



v a - v b 

In the case of one planet in circular orbit and 
homogeneously distributed observations, (Wb) = 
K 2 /2 (see the Appendix and Cumming 2004). 

The only random quantity in equation (|8|) is 
Qb- It is known from linear regression theory (see 
Kurth, 1967) that Qb/cr 2 is a xl variate with v = 
Vb. Therefore, the statistic z\ is an inverse chi- 
square distribution, and 

Prob(z > Zd) — Prob (z\ > z\d) — Prob ( — > z\d 

\ Xii 

(11) 

From these expressions, we can obtain the mean 
and variance of z as 

A 



E(z) 



D\z) 



Vb 



Vb 



V a - V b 



(12) 



2A 2 



K-2)> b -4) 



2 Throughout the paper the symbol over a given quantity 
indicates its best-fit estimation. 
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Fig. 3. — Histogram of the best-fit values z result- 
ing from 774 random data sets constructed for a 
single planet (mi = lMj up ) in a circular orbit with 
orbital period of 100 days. For each sample the ini- 
tial mean longitude was chosen randomly, and the 
corresponding RV curve was sampled with iV = 60 
and K/a = 3. The solid curve shows the inverse 
chi-square distribution using equation (jTTJ) , while 
the dashed line corresponds to the Gaussian dis- 
tribution with same mean and variance. All plots 
subtend a unit area. 

In the approximation for large N — Mb, the stan- 
dard deviation of z is approximately given by 



V2N 



(M b - l)y/N-M b 



(W fc 2 > 



(13) 



These results may be compared with those pub- 
lished by Cumming (2004). If we adopt the ap- 
proximations (W 2 ) = K 2 /2 and a e — a, the re- 
sults for (z) are not very different from those given 
above and tend to be equal for very large N. How- 
ever, Cumming's results for cr| are significantly 
smaller than those given by equation (fT2"]) (e.g., 
9.8 instead of 13.7 in the example considered in 
Fig. E|) . The corresponding Gaussian is then much 
more narrow than than shown in Fig. [3] and does 
not reproduce satisfactorily the distribution of the 
values of z. 

For a given signal amplitude K, the distribu- 
tion of z gives the dispersion of possible powers 
obtainable from different data sets with the same 
N and signal-to-noise ratio. Thus, for given de- 
tection threshold Zd (or z%d), we can look for the 



value of the parameters corresponding to a given 
probability of detection Pdotoct- This is done solv- 
ing the equation Prob \Jg > Zidj = F d ctcct with 
respect to z\d- Adopting a numerical value for 
^detect (e.g., 99%), equations (fTT )) -([T2 ]) yield the 
necessary value of A. Finally, from the approx- 
imation (yV 2 )/of = K 2 /(2a 2 ) we can obtain a 
relation between K/a e and N to satisfy this crite- 
rion. It is worth mentioning that the results thus 
obtained are very robust and do not depend criti- 
cally on the chosen distribution of the best-fit val- 



3.2. Detectability of Two Resonant Plan- 
ets 

To estimate the detectability limit for two plan- 
ets in the 2/1 MMR, we assume that the body gen- 
erating the largest amplitude in the RV signal is 
already known, and wish to calculate under which 
conditions the signal of the second body also satis- 
fies the given detectability limit. A look at the dif- 
ferent plots of Fig. [2] shows that there is no unique 
separation of the signals. When the mass ratio 
1712/ mi is large, the RV amplitude of the outer 
planet is also large (implying K 2 > Ki), and the 
detectability of the resonant pair reduces to the 
discernability of the inner mass. In the opposite 
case where K\ > K2, it is the inner planet that 
dominates the RV curve and the existence of the 
outer body must be deduced from the difference 
between successive maxima. 

Extending the equations in Section [2~T1 to ellip- 
tic orbits, we can write the ratio of amplitudes of 
both planets as 




(14) 



showing that the value m^/mx ~ v2 deduced ear- 
lier still provides a qualitative critical value for 
K\ ~ K2 when both planets have similar eccen- 
tricities. 

We begin discussing the case of a more mas- 
sive inner planet. Here K\ > K 2 , and we wish 
to estimate the detectability limit of finding ni2 
assuming that mi is known. We can employ ex- 
actly the same criterion constructed for one-planet 
systems, where the number of free parameters are 
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now Mb = 11 and M a = 6, and rewrite 134] 



vb /||W 2 || + (iV-l)^ : 



Va - Vb 



(15) 



If we assume a' e ~ cr (see the Appendix), the re- 
maining procedure is analogous to that deduced 
in Section 13.11 It is then possible to estimate a 
limiting value of Kij 'cr, as function of the number 
of points N necessary for the resulting peak to 
surpass a pre-established FAP value with a prob- 
ability equal to Pdetect- Moreover, since 



K2 

a 



K2 



Ki 
a 



(16) 



we can also address this problem saying that the 
detectability condition depends on N and on the 
ratios K\ju and K2/K1. Applying equation (fl"4"|) 
we can transform the condition in K2/K1 to one 
in m-i/rn\. Since Ki < K\, the question reduces 
to finding what values of ro2/mi (for given eccen- 
tricties e») reduce the value Kijo to undetectable 
magnitudes. 

The ratio K\ / cr depends on the standard devia- 
tion of the data, the stellar mass, as well as on the 
orbital parameters and mass of the inner planet. 
For example, larger values of mi will increase the 
amplitude of K\, thus allowing the detection of 
smaller values of mijmx for a fixed a. Let us call 
Ki the RV amplitude generated by a Jovian size 
planet (mi = Mj up ) at a\ = 1 AU and eccentric- 
ity e\. Since K\ oc nia\(mi / 'too) , for any other 
mass and semimajor axis, the ratio K\/a may be 
written as 

Ki _ K 10 
a cr* ' 

where the "scaled" standard deviation cr* 
lated to the nominal value a by 



(17) 



is rc- 



/ mi 



\M Jup J VlAU 



«i 



-1/2 



— ] • (18) 



Finally, from equations (fl4|) . (|T6|) and (JXTJ) we can 

then write 




m 2 



1 — e| TO-i 



Kip 
a" 



This equation allows us to estimate the minimum 
value of TO2/TO1 for detectability of the outer res- 
onant planet as function of N and cr*, for given 
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Fig. 4. — Minimum value of mass ratio to 2 / TOx 
necessary for detectability of the smaller outer res- 
onant planet (with Pdetect = 0.99), as a function 
of number ./V of data points and for four values 
of the scaled standard deviation cr*. Both plan- 
etary eccentricities are chosen as e, = 0.1. The 
FAP is taken equal to T = 10~ 4 . The continu- 
ous lines correspond to predictions with the in- 
verse chi-square distribution pT|) . while dashed 
curves are obtained with the Gaussian approxi- 
mation with equal mean and variance. 

orbital eccentricities e,. If we assume that the in- 
ner planet corresponds to a Jovian mass at 1 AU, 
cr* is equal to the nominal standard deviation of 
the RV data a. For other inner masses, cr can be 
calculated through (fT8]) . 

Fig.|4]shows the detectability limit for m^l mi < 
1 for several values of a* and as a function of the 
number N of data points. Each curve gives the 
minimum value of the mass ratio to 2 /toi neces- 
sary to achieve a FAP of T = 10" 4 . The stellar 
mass was too = 1M Q and both eccentricities were 
chosen equal to 0.1. Continuous lines were con- 
structed using the inverse chi-square distribution 
pip while the dashed curves denote approximate 
values obtained with a Gaussian distribution with 
same mean and variance. Qualitatively both mod- 
els yield to similar results, although expression 
(fTTjl leads to more restrictive detectability limits. 

For cr* = 15 m s _1 , practically no resonant sys- 
tem is detectable with mass ratio to 2 /toi lower 
than ~ 0.6, even for a data set containing N — 200 
RV observations. For N < 60 practically no sys- 
tem with this cr* is detectable for any mass ra- 
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Fig. 5. — Detectability of two planets in the 2/1 
MMR. Open circles show results of Monte Carlo 
simulations from synthetics data sets. The plot 
show the FAP for random values of the mass ratio 
m-2/mi. The total number of points was N — 100 
distributed randomly over four orbital periods of 
the outer body. The broad continuous curves show 
prediction values from the model with Pdetect = 
0.99. 

tio. The picture improves with lower values of the 
scaled standard deviation, until for a* = 4 m s _1 it 
is possible to detect systems with mass ratio close 
to one-tenth provided the data set is sufficiently 
large. 

It is important to note that the value of a to 
consider in these calculations is not just the one 
deduced from the observational techniques, but 
the total standard deviation including stellar jitter 
and the possible incompleteness of the two-planet 
model. Typical values of stellar jitter (Wright 
2005) are of the order of 4 — 7 m s _1 , and instru- 
mental uncertainties are also of similar magnitude. 
The square sum of both implies that typical val- 
ues of <7 for these systems should be of the order 
of 6 — 10 m s _1 . A good qualitative estimate of 



a may be given by the value of the rms resulting 
from the orbital fit (Shen & Turner 2008) which 
leads to a ~ 7 — 10 for most multiple-planet sys- 
tems. 

To check the precision of the limits obtained 
with this model, we performed a series of Monte 
Carlo simulations with synthetic data sets, each 
with N = 100 data points covering a total of four 
orbital periods of the outer resonant body and 
with mass ratios in the interval mijvnx € [0.1, 1.0]. 
The inner planet was fixed to the same mass and 
orbital parameters as used in Fig.|H(mi = IMjup, 
a\ = 1 AU). One series of simulations was done 
adopting a* = 4 m s , while a second run con- 
sidered it = 10* m s — 1 . For each data set we 
determined the best two-planet fit and calculated 
the value of the FAP corresponding to the sec- 
ond planet using the same procedure as for the 
one-planet system. Results are shown in Fig. 
where the open circles show the relationship be- 
tween mass ratio and the respective value of FAP 
for all the data sets. Values of the ordinate lower 
than 10 -13 were equated to this lower limit. Fi- 
nally, the solid curves present the analytical de- 
tectability limits obtained from the model. 

These plots must be read in the following way. 
For a given value of the mass ratio 1112/011, all the 
synthetic data sets give orbital fits with FAP lower 
than the broad continuous curve. Thus, in order 
to guarantee detectability, the FAP corresponding 
to a given mass ratio must be above that thresh- 
old. The agreement between the analytical curve 
and the numerical simulations is very good, indi- 
cating that the simplified model described above 
is more than adequate to predict the detectability 
of multiple planet system. 

The same calculations can be done for systems 
with K2/K1 > 1. In such a case the detectabil- 
ity criterion follows the same route, but now the 
outer planet is responsible for the dominant peak, 
while the signal that defines the detectability of 
the resonant system is K\ . Consequently, equation 
(|16p must now be inverted and will yield the mini- 
mum necessary value of K1/K2 such that the value 
K\Ju is still detectable. The dependence with the 
individual masses and semimajor axis can again be 
overcome defining the "scaled" standard deviation 
a* as in equation (p~8|) . Note that we have kept the 
mass and semimajor axis of the inner planet as the 
scaling parameters. This is a matter of choice but 
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Fig. 6. — Same as Fig. [4J but including mass ra- 
tios m^/mi > 1. Only results constructed with 
the inverse chi-square distribution (llip are shown. 
For a given value of a* , detectability is guaranteed 
only inside the two corresponding curves. 

helps maintain a certain homogeneity, regardless 
of the mass ratio under consideration. 

Fig. [6] shows the detectability limit combining 
results for all mass ratios. For a given a*, all de- 
tectable planetary systems lie inside a pair of cor- 
responding curves. Note that the intersection be- 
tween both sets of curves occurs for 7712/mi ~ s/2, 
as expected from the condition K2 = K\ . Finally, 
there is a certain variation with the orbital eccen- 
tricities stemming from the relation between sig- 
nal amplitude and planetary mass (see equation 
(fT4| ). However, qualitatively the picture remains 
the same. Except for very small standard devia- 
tions, it appears very difficult to detect resonant 
planetary systems with mass ratios m 2 /mi larger 
than ~ 4 or lower than ~ 0.3. It is important to 
stress that this does not mean that it is impossible 
to detect systems outside these borders; but this 
detectability will not be guaranteed for any data 
set and will depend on the specific time distribu- 
tion and/or errors of individual data points. 

The previous analysis assumes that m\ = 
Aij up , a x = 1 AU and m = M Q , for which 
a* = a. Other values however, will change the pic- 
ture. Among the known resonant population, per- 
haps the most extreme case is the GJ876 system 
(1712/1711 ~ 3.1), where mi = 0.6Afj up , a\ = 0.13 
AU and m = 0.32 Af Q (Butler et al. 2006). Ap- 
plying equation (fTSl) to the results shown in Fig. [6j 
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Fig. 7. — Detectability limits of the mass ratio 
m2/i7ii (with Pdetect = 0.99 and N = 100), as a 
function of the inner mass mi, for four values of 
the standard deviation a . Both planetary eccen- 
tricities are chosen as = 0.1. The FAP is taken 
equal to T = 10~ 4 . As before, for a given value 
of a, detectability is guaranteed only inside the 
two corresponding curves. Above the broad black 
curve 77i2 > 15A/j up , a value usually associated 
with the brown dwarf limit. 

we find that for N = 100 the resonant system can 
be identified even for standard deviations in the 
RV data of the order of a ~ 40 m s _1 (correspond- 
ing to er* ~ 8 m s _1 ). 

Fig. [7] shows how the detectability limit of 
7712/777,1 varies with the magnitude of mi, now as- 
suming a fixed length of the data set N — 100. 
Both ai and itiq remain at their original values. 
Once again, a resonant system is detectable if the 
mass ratio 777,2/7774 is located inside two curves 
of equal color. We can see that the detectabil- 
ity increases significantly with the mi, allowing a 
larger range of mass ratios than previously shown. 
However, the possibility of very massive resonant 
planets must be considered with care since the 
augmented mutual perturbations may compromise 
the dynamical stability of the system. 

The difficulty in detecting resonant planets 
with significantly different masses is in accordance 
with currently known systems. Recall that Fig. [T] 
shows the distribution of mass ratios against or- 
bital period ratio for all consecutive planetary 
pairs in multi-planetary systems. All planetary 
systems near the 2/1 MMR have mass ratios 
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Fig. 8. — Distribution of ratio of mean motions 
as function of planetary mass ratio m,2 jm\ , from 
the Monte Carlo simulations used in Fig. [5l Only 
those systems with T < 10 -4 were used. 



77i2 /mi grouped near unity, in accordance to the 
analytical findings of the present model. The in- 
teresting conclusion is that it is possible that reso- 
nant systems may exist for other mass ratios, but 
are presently undetectable. 

4. Error Estimation in the Fitted Param- 
eters 

Even if the two-planet resonant system is ac- 
tually detected, there is no guarantee that the 
masses and orbital parameters can be estimated 
with any accuracy. As an example, Fig. [H] shows 
the distribution of ratios of mean motion for those 
synthetic data sets in Fig. [5] that satisfied the de- 
tectability criterion. A significant dispersion in 
the mean motions is noted around the nominal 
value ri\jni — 2, with an appreciable percent- 
age of the systems which would be qualified as 
near-resonant but not locked in MMR. This effect 
presents an additional difficulty in detecting reso- 
nant planets for mass ratios m^/mi different than 
unity. 

These results were obtained with synthetic data 
sets covering four orbital periods of the outer 
planet. Although the detectability criterion is in- 
dependent of the total observational time interval 
(as long as we can guarantee a good coverage of 
the phases), the precision of the orbital elements 
is expected to improve with longer time intervals. 
Nevertheless, even if correctly identified as being 
locked in resonance, the planetary orbital elements 
may still contain significant errors and affect the 
deduced dynamics for the fitted systems. To quan- 
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Fig. 9. — Results of best fits of 1000 synthetic data 
sets corresponding to a nominal solution given by 
the first data set in Table [fl with N = 200 data 
points and constant standard deviation of a = 7 m 
s -1 . The large filled gray circles indicate nominal 
solution. In the eccentricities plane the gray curve 
marks the family of zero-amplitude ACR solutions 
for a mass ratio 7712/7711 =3/1. 



tify this effect, will analyze three test cases, one 
with mass ratio larger than one, one with equal 
masses and one with mass ratio lower than unity. 

4.1. More Massive Outer Planet 

We begin analyzing the case where the outer 
planet is more massive than its inner companion 
(i.e., 7712/mi > 1). We have chosen the first set 
of initial conditions shown in Table [TJ character- 
ized by a mass ratio of 7712/7711 = 3/1 and orbital 
elements leading to a small amplitude (~ 5°) os- 
cillation around a (0, 0)-type ACR. We generated 
1000 synthetic RV data sets, each with N = 200 
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data points distributed randomly over four orbital 
periods of the outer planet (i.e., four years). The 
values of each individual RV were taken randomly 
following a Gaussian distribution around the exact 
value and a variance given by a constant standard 
deviation for u = 7 m s~ 4 (i.e., a* = 5.5 m s^ 1 ). 
These values guarantee a detectability of the two 
planets and places the deduced mean-motion ratio 
within the commensurability region. 

For each data set we first calculated a single- 
planet orbital fit starting with a genetic algorithm 
(Charbonncau 1995) without any pre-established 
initial solution. If the resulting fit had an as- 
sociated FAP below 1CP 4 it was accepted and a 
two-planet fit was performed. The single-planet 
solution was used as a first guess and the code 
searched for the two-planet solution with the min- 
imum rms. Once again this was accepted only if 
the FAP of the two-planet fit with respect to the 
single planet solution was below 1CP 4 . If accepted, 
we then numerically integrated the configuration 
over a time span of 10 4 orbital periods and esti- 
mated the amplitude of oscillation of the resonant 
angles 9i = 2X2 — \\—vj{ and the difference in lon- 
gitudes of pericenter Aw = wi — w\ . The critical 
angle 9\ corresponds to an interior resonance (see 
Michtchenko et al. 2008a), as found in symmetric 
ACR, while 82 is the librating angle for exterior 
resonances (see Michtchenko et al. 2008b) and is 
thus applicable to asymmetric ACR. 

Results are shown in Fig. [51 where we present 
the distribution in several planes. In each frame 
the large filled gray circle marks the nominal solu- 
tion and the open black circles the best fits of each 
synthetic data set. The gray curve in the eccen- 
tricities plane shows the (0, 0)-family of ACR for a 
mass ratio 7712/7711 equal to 3/1. As expected from 
the choice of N and a, all fictitious two-planet sys- 
tems were effectively detected by the fitting pro- 
cess and also identified as being locked in the 2/1 
MMR, and both 9\ and Aw oscillated around the 
equilibrium solution. However, on 2% of the cases 
Aw circulated taking all values between and 2n, 
while 6\ librated. This type of motion is usually 
referred to as a #i-libration and is considered dif- 
ferent from an ACR. However, as shown recently 
by Michtchenko et al. (2008a), even within an 
ACR the libration of Aw is purely kinematical 
and not associated with any separatrix crossing as 
long as the eccentricities are not too high. 
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Fig. 10. — Distribution of calculated eccentricities 
and amplitudes of oscillation of the resonant angle 
9 1 and Aw for the case 7712/7711 =3/1. Nominal 
values are indicated by the vertical dashed lines. 
The color codes are used to identify different sizes 
N of the data sets. 

The two top graphs in Fig. [9] show the distribu- 
tion of planetary masses. From the left-hand plot 
we can see that the larger mass (i.e., larger K) has 
a fairly small dispersion of values, although mi is 
not very well determined leading to a large disper- 
sion in the calculated mass ratio mijm\ (right- 
hand graph) between 2.5 and 3.5. Even so, the 
semimajor axis and ratio of mean motion are very 
well estimated with little variation around the ex- 
act values (mid left-hand plot). 

The distribution of solutions in the eccentric- 
ities plane shows an interesting correlation with 
the family of zero amplitude ACR (shown in gray). 
This is an unexpected result since there is no ev- 
ident relationship between the fitting process and 
the topology of the MMR. However, it does in- 
dicate that even imprecise determination of the 
eccentricities will not necessarily lead to large 
displacement from the ACR family and, conse- 
quently, from stable solutions. Recall that the 
successive orbital determinations of the resonant 
GJ876 planets (mass ratio close to 7712/7711 = 3/1) 
show a similar behavior. Although in the past 
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eight years the eccentricities have varied signifi- 
cantly, all best fit solutions have been found close 
to a stable ACR. 

Finally, the two bottom frames show the vari- 
ation in the angular variables. The graph on the 
left gives the distribution of the calculated values 
of 61 and Azu. More interestingly, the right-hand 
frame shows the distribution of the amplitudes of 
oscillation around the ACR. Note a very impor- 
tant dispersion in the dynamical behavior of the 
calculated fits. As with the eccentricity, the am- 
plitude has a lower bound at zero which causes a 
systematic effect leading to the fact that the es- 
timated amplitude of libration is always overesti- 
mated with respect to the true value. 

This figure shows that, even under favorable 
circumstances given by a large data set and rea- 
sonable observational uncertainties, the errors in 
the orbital fit are still significant, especially in the 
eccentricities and angular variables. Even if the 
ACR is still identified as such, the amplitudes of 
oscillation are greatly increased beyond their nom- 
inal values. 

Fig. [TO] shows the distribution of fitted eccen- 
tricities and amplitudes of oscillation of the angu- 
lar variables for N — 200 (black histograms) and 
N = 100 (red histograms). The nominal values 
corresponding to the original system are marked 
with the vertical dashed lines. In all cases the data 
points were distributed randomly over four orbital 
periods of the outer planet and a = 7 m s _1 . 

The two top frames correspond to the eccentric- 
ities; ei on the left and ei on the right. The distri- 
bution of eccentricities of the inner planet is fairly 
symmetric with respect to the true value (~ 0.42), 
and is reminiscent of the results shown for sin- 
gle planet systems (Shen & Turner 2008). Al- 
though the dispersion increases slightly for smaller 
N, there appears to be little systematic error and 
the peak of the histogram is close to the vertical 
dashed line. The same, however, is less evident 
in the case of the eccentricity of the outer planet, 
where the distribution is less symmetric and there 
is a bias towards larger values of ei . 

The two lower plots give the distribution of am- 
plitudes of oscillation of the resonant angle 9\ (left- 
hand graph) and Atn (right-hand graph). Once 
again the nominal values are shown in the vertical 
dashed lines and correspond to a small amplitude 




Fig. 11. — Same as Fig. [TO] for mass ratio 
tail m\ = 1/1- However, now color codes are used 
to identify different values of mi, maintaining a 
constant length of the data set N = 100. The 
black histograms were obtained with mi = lMj up 
while those in red were calculated with mi = 
2M Jup . 

solution. As shown in Fig. [5] there is an important 
increase in the amplitudes leading to a systematic 
bias towards solutions with larger amplitudes of 
oscillation. Surprisingly there is little change for 
smaller data sets. In the case of 0\ there is practi- 
cally no calculated solution with amplitude larger 
than 90°. The same, however, does not hold for 
Aw. Even for larger data sets, some solutions 
have amplitude equal to 180 degrees correspond- 
ing to circulations of the difference in longitudes 
of pericenter. 

4.2. Equal Mass Planets 

The previous analysis can be repeated for a 
mass ratio of unity (7712/mi = 1) and choosing 
initial conditions close to a stable ACR. Orbital 
elements are once again indicated in Table [JJ This 
particular ACR was chosen such that the total an- 
gular momentum of the system was approximately 
the same as in the previous case; mainly this im- 
plies that the sum of the square of the eccentrici- 
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mijm\ 




Non resonant (%) 


fli-Librator (%) 


ACR (%) 


100 




3/1 


N = 200 





2 


98 


75 






N = 100 





12 


88 




1/1 


mi = 1 


9 


21 


70 


50 






mi = 2 


2 


2 


96 






1/3 


mi = 1 


68 


6 


26 


25 






mi = 3 


23 


1 


76 







Table 2: Percentage of Different Dynamical Out- 
comes of the Best Fits for Different Mass Ratios 
7712/7711. For symmetric ACR the resonant angle 
is taken as 9\, while #2 is adopted for asymmetric 
solutions (see Michtchenko et al. 2008ab). The 
values of mi are given in units of Mj up . 

ties is of the same magnitude. 

Fig. [TT] shows histograms with the distribution 
of eccentricities (top), amplitude of libration of 9\ 
(bottom left) and amplitude of oscillation of Avj 
(bottom right). Contrary to Fig. QUI the number 
of data points has been kept constant at N = 100 
and the color code now identifies different values 
for the mass of the inner planet. The black his- 
tograms were drawn using mi = Mj up , while for 
the plots in red we adopted mi = 2Mj up . This 
latter value generates a RV amplitude in the com- 
plete signal similar to the case discussed in Section 
4.1, thus allowing an easier comparison of the dis- 
persion of the orbital solutions. 

As before, the eccentricity of the inner planet is 
fairly well established. For both values of mi the 
distribution is centered near the nominal value, al- 
though we note a smaller dispersion in the results 
for larger masses (i.e., larger amplitudes of the RV 
signal). The eccentricity of the outer planet, how- 
ever, shows a more complex distribution, with a 
notorious systematic shift in the averaged value. 
Also as before, practically all the orbital fits lead 
to large amplitudes of libration, especially of the 
difference of longitudes of pericenter. The mean 
value of this last angle lies near 90°. Although 
these bias decrease significantly with larger values 
of mi, they are still much larger than those pre- 
sented in Fig. [lOl 

4.3. Less Massive Outer Planet 

Finally, we analyze the case where 7712/777.1 < 1. 
Initial conditions (Table [TJ correspond to a small 
amplitude oscillation around an asymmetric ACR 




Fig. 12. — Same as Fig. [TT] for mass ratio 
7712/7711 = 1/3. Black histograms were obtained 



with mi 
3M Jup . 



lAfjup while those in red with mi 



with 7712/7711 = 1/3. Once again the histograms in 
red were obtained with increased planetary masses 
(mi = 3-Mjup) which yield a magnitude in the 
complete RV signal similar to the initial conditions 
for 7712/7711 = 3/1. 

Although both planets were detected in all the 
data sets, for mi = Mj up only 26% of the resulting 
orbits reproduced the ACR. In 68% of the cases 
the resulting configuration was non-resonant and 
neither 62 nor Azu librated around the equilibrium 
values given by the ACR, but instead circulated 
taking all values between zero and 27r. Given the 
high eccentricity of the outer planet, these orbits 
inevitably led to unstable motion and a disruption 
of the system in short timescales. In the remaining 
6% of the cases the best fits yielded # 2 -hbrations 
but a circulation of the difference in longitudes of 
pericenter. 

These numbers improve significantly for mi = 
3Mj up , where disruptions occurred only in 23% 
of the cases and most of the rest corresponded to 
stable ACR. Even so, this implies that for scale 
standard deviations as low as a* = 2 m s _1 , a 
quarter of the synthetic data sets failed to iden- 
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tify the resonant motion and led to dynamically 
unstable fits. 

The distribution of the orbital elements of the 
stable ACR solutions is shown in Figure [T3J Since 
we have eliminated the non resonant cases, the 
histograms appear similar to the other mass ra- 
tios, with a fair estimation of the eccentricity 
of the inner planet and a much larger disper- 
sion in ei. For m\ = -Mj up we also note indi- 
cations of a bimodal distribution for the eccen- 
tricity of the outer planet, with one peak cen- 
tered approximately around e 2 = 0.45 and a sec- 
ondary peak near 0.6. This bi-modality is related 
to the bifurcation of the family of asymmetric 
ACR for this mass ratio and high eccentricities 
(see Michtchenko et al. 2008b). However, the his- 
togram in red shows a more localized distribution 
of e 2 . 

As with the other mass ratios, the distribu- 
tions of amplitudes of libration/oscillation also 
show a marked bias towards larger values. Al- 
though this decreases for larger values of mi , even 
for mi = 3Mj up practically none of the best fits 
are able to reproduced the small-amplitude oscil- 
lation around the asymmetric ACR. Finally, in 
the black histogram the distribution of the am- 
plitude of libration of 62 shows a marked peak 
around 60 degrees, which roughly corresponds to 
the maximum amplitude permitted by the asym- 
metric ACR. Larger amplitudes place the system 
outside the resonant domain, which may explain 
the large proportion of non-resonant outcomes in 
the orbital fits. 

4.4. Application to Real Exosystems 

Table [2] summarizes the possible dynamical out- 
comes of the orbital fits for the three mass ratios. 
Although all data sets were constructed with the 
same standard deviation a = 7 m s _1 , there is a 
marked difference in the results. For m-ijm\ =3/1 
most fits correctly identified the ACR and all led 
to stable motion of both planets in timescales of 
the order of 10 4 orbital periods. The inverse case, 
Tn&fmi = 1/3 is significantly less precise, and is 
probably due to the existence of several domains 
of stable motion, each associated with a different 
family of ACR (Michtchenko et al. 2008b). Each 
domain is separated by a chaotic layer, and are 
less robust than those found for mass ratios larger 
than unity. 




Fig. 13. — Distribution of planetary eccentrici- 
ties from best fits of synthetic data sets, compared 
with the families of ACR. The color codes identify 
different mass ratios. Black : m^/mx — 3/1, red : 
m^/mi = 1/1 and blue : mz/mx = 1/3. (a): All 
inner planetary masses equal to lMj up . (b): Aug- 
mented inner planetary masses: mi = 2Mj up for 
m,2/mi = 1/1 and mi = 3Mj up for rri2/mi = 1/3. 

It thus appears that resonant planets with mass 
ratios lower than unity are not only more diffi- 
cult to detect, but even when detected the orbital 
fits usually lead to wrong period ratios and/or 
dynamically unstable solutions. Although these 
results are valid for the 2/1 resonance and the 
55Cnc-c and 55Cnc-d planets are located in the 
vicinity of the 3/1 commensurability (mass ratio 
7712/ 'mi ~ 1/5), it seems plausible to expect simi- 
lar results, particularly since asymmetric ACR are 
involved in both cases. This could be a possi- 
ble explanation of why successive orbital fits of 
the 55Cnc planets alternatively point towards res- 
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onant or non-resonant motion, while the same un- 
certainty is not present for exosystems with other 
mass ratios. Another example is given by the 
HD37124 system with two planets in the vicinity 
of the 2/1 MMR. Although the best fit indicates 
mijm\ ~ 1.2 and an asymmetric ACR (Baluev 
2008), other solutions exist in which the bodies 
are outside this commensurability. 

Fig. [T3] compares the distribution of eccentrici- 
ties found for mijm\ =3/1 (black open circles), 
m%jm\ = 1/1 (red) and rr^/mi = 1/3 (blue). 
Only those orbital fits leading to stable ACR were 
plotted. The colored curves correspond to the 
zero-amplitude ACR for each mass ratio. The top 
plot was drawn using mi = Mj up in all three cases. 
For the bottom plot we used the augmented values 
of mi given in the caption. 

For m,2 /mi = 1/3 we note that most of the so- 
lutions lie relatively close to the family of stable 
ACR, although it is important to recall that these 
correspond to less than half of the total data sets. 
A similar trend is observed for m<i/m\ — 3/1, ex- 
cept that all orbital fits were now plotted. How- 
ever, in the case of equal masses 1712/mi = 1/1, 
there is no apparent correlation between the ec- 
centricity dispersion and the ACR families is ab- 
sent in this case. Thus, it appears that, for large 
mass ratios, there is a much larger probability that 
even uncertain orbital fits will remain close to sta- 
ble ACR configurations. For mass ratios close 
to unity, this is in no way guaranteed, and even 
small standard deviations may lead to larger am- 
plitude of libration or unstable motion. Even with 
increases planetary masses and, thus, better or- 
bital fits, the same correlation is still found. In 
fact, the large signal-to-noise ratio of the RV data 
maintains a similar dispersion in e%, although the 
estimated value of ei is more precise. 

5. Conclusions 

We have analyzed the detectability limit of two 
planets in the vicinity of a 2/1 MMR, as a func- 
tion of the number N of data points and the stan- 
dard deviation a of the RV values. For values of 
a and N comparable to real exosystems we found 
an important bias toward cases with mi/mx ~ 1- 
Planets with mass ratios much larger or smaller 
than unity are much more difficult to separate 
from Doppler data. This may not only help ex- 



plain the lack of resonant systems with these mass 
ratios but also may indicate that the number of 
resonant systems in existence is actually signifi- 
cantly larger than currently believed. 

We also analyzed the dispersion in the masses 
and orbital elements of the detected systems. Al- 
though for large mass ratios most synthetic data 
sets effectively detect a resonant motion, for low 
values of 1712/1111 a large proportion of the or- 
bital fits fail to detect the resonance lock and er- 
roneously indicate a non-resonant configuration. 
Even within the subset of resonant orbits, we note 
a large dispersion in the estimated values of the 
eccentricity of the outer planet, and the median of 
the distribution is always higher than the nomi- 
nal value of This is in accordance with sev- 
eral real exosystems in the 2/1 MMR, particu- 
larly HD82943. Moreover, nominal solutions cor- 
responding to small-amplitude librations around 
ACR are usually accompanied by an important 
bias toward incorrect large amplitudes that may 
even compromise their dynamical stability on long 
timescales. 

Even if this is not the case, it presents a pos- 
sible answer to the question of why, with the ex- 
ception of the resonant bodies in GJ876, all other 
resonant planets have orbital fits consistent with 
large amplitude librations. Although recent works 
by Sandor, et al (2007) and Crida et al. (2008) 
propose that the observed large amplitudes of li- 
bration may have been caused by past dynamical 
processes (e.g., planetary scattering or disk disper- 
sal), the present results indicate that this observa- 
tional characteristic may not be real, but a simple 
consequence of dealing with small data sets with 
large uncertainties in the RV values. 

For multiplanetary systems it is important to 
stress that the absolute numerical precision of the 
orbital elements is not, in itself, the only datum 
to consider, as much as the diversity in dynamical 
behaviors within the uncertainty region. In other 
words, if the two planets are non resonant with 
low eccentricities, even significant changes in the 
orbital elements will not necessarily lead to differ- 
ent regimes of motion. Conversely, if the planets 
are locked in a MMR and (preferably) have mod- 
erate to high eccentricities, even a small change 
in the orbital parameters may mean the difference 
between stable and unstable motion. 
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Appendix. Data fitting to a model with or- 
thonormal functions. 

We assume that the adopted model W(t) is 
given by a linear combination of M independent 
continuous functions H^t) (// = 0, 1, • • • , M — 1). 
This is, for instance, the case of the very simple 
model of one planet in circular orbit, for which the 
true and mean anomalies are equal / — £ (for a 
given fixed period) . Although the general (elliptic 



This 2-column preprint was prepared with the AAS IATjrjX 
macros v5.2. 
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orbit) model is no longer linear, we can linearize 
it in the neighborhood of the true values, which 
are then assumed to be determined by successive 
iterations. We assume that the functions H^(t) 
can be orthonormalized to generate a new set of 
M functions such that 

(VM = <W ( 2 °) 

where 5 MiP is the Kronecker symbol and {h^.hp) 
represents the inner product 

N 

{h IM h p ) = J2 h »( t i) h p( t i)- ( 21 ) 
1 

If the observations do not have the same variance, 
weights may be introduced in this definition to 
take it into account. Hence 

M-l M-l 

w(u) = /w**) = E a ^ ( 22 ) 



(see Kurth, 1967). The construction of the or- 
thonormal functions is a cumbersome step. In the 
case of the search of periods of unequally spaced 
observations, they were first introduced by Ferraz- 
Mcllo (1981). In this application, it is not neces- 
sary to actually construct the orthonormal func- 
tions h(ti). It is only necessary to know that such 
operation is possible and that these functions al- 
low much easier derivations as they behave just as 
orthogonal unit vectors. In the same way, we do 
not need here the actual definition of the a's (as 
functions of the /3's). 

One important point often overlooked is that 
the set of function HJt), for reasons of complete- 
ness, must include the function H (t) = 1. In 
Scargle's (1982) periodogram, for instance, this 
function is not included. However, the assump- 
tion that the given data have zero average (i.e., 
J2i y{U) = 0) docs not guarantee that the data 
modeled at the given times will also have zero 
average. In general, if W(U) are the radial ve- 
locities computed at the times U using the pro- 
posed model, we have ^ 0. Conse- 
quently an incomplete basis, in which H (t) = 1 
is missing, cannot allow the function that model 
the data to be represented. When just the pe- 
riod of such a function is looked for (the usual 
question that a periodogram analysis is assumed 
to answer), a basis formed by sin(wt) and cos(uit) 



may be enough, but for the full problem of best 
fitting one sinusoid to the data, it is not. For this 
reason, in many applications, as those concern- 
ing "detectability" using periodograms, the com- 
plete model given by Fcrraz-Mcllo (1981) (or some 
equivalent ones published later; see references in 
Cumming et al. 1999) must be used instead of the 
incomplete ones. 

The following results are elementary: 

• The solution Q = Q m i n of the best fit prob- 
lem is given by 

3„ = (V,fc„)> ( 23 ) 

which is obtained from dQ/da^ = 0. Note 
that ho = 1/N and, therefore, So = (V) ; 

• In the best fit, the square sum of the resid- 
uals is given by 

Q = \\V\\-\\W\\ (24) 

where 1 1 V| | = (V, V) and | \W\ | is the spectral 
power 

M-l 

\\W\\ = (W,W)=J2^ (25) 
o 

• If we assume that the data are Vi — W(U) + 
6i where W is a given model function and 

a normal variate N(0, a 2 ) uncorrelated with 
W (i.e., (W,e) = 0), there follows ||V|| = 
1 1 W| | + | |e|| and, therefore, 

11^11 = HW|| + ||e||-Q; (26) 

• If the observations are perfectly distributed, 

||W||=£n? = AT<W 2 >= lim - / W(t) 

(27) 

In the case of one planet in circular orbit, 
this operation gives the value NK 2 /2 found 
in the quoted papers. However, even in this 
more simple case, if the observations are not 
well distributed, the non-orthonormality of 
the functions H^(t) leads to a different re- 
sult. 
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Note that when the given data averages to zero: 
(V) = 0. This often occurs because it is usual in 
periodogram studies to refer the given data to its 
average. In such case So = and ||V|| = (V,V) is 
the variance of the given data (Vi) times (N — 1). 

Two particular cases are important in the ap- 
plication to exoplanetary detection. 

1. One-Planet Case: 



Here model a is a constant (i.e., W a — fy). 
Then we can simply write 1 1 W a \\ = a$ = 
(V) 2 . If, in addition, the given data are such 
that (V) = 0, then the floating mean power z 
(see equation (QJ) corresponding to the best 
fit solution is 



g. _ N — M b \\W b \\ 
M b -I Q b 

Introducing (j2€3[) we obtain 



AT -M 6 /||W b || + ||e|| 



M b -1 V 



(28) 



(29) 



Considering that ||e|| = e 2 is an estima- 
tion of (N — 1) a 2 , where a e is a characteris- 
tic standard deviation of the RV data, then 
finally 



N - M b / \\W b \\ + (N - 1)<j 2 
M b -I \ 



Qb 



l 



(30) 



where recall that W b is the signal generated 
by the single planet. 

2. Two-Planet Case: 

We now assume that model a is the one- 
planet model and b the two-planet model. 
As before, the synthetic sample is given by 
V{ = yV b (U) + £i where is a Gaussian noise. 
We assume that W b (ti) is the sum of two 
one-planet models W\{ti) and W2{ti) 1 where 
the index are chosen such that the signal 
of the first planet is larger (i.e., ||Wi|| > 
1 1 VWa ID- We can then write 

IM| = ||Wi|| + ||W 2 || + (Wi,W 2 )+||e|| (31) 

where (VV^Wy is the covariance between 
the signals generated by each individual 



planet. Introducing this expression into 
we can write 



-Wa\\ 

(32) 



Qa = ||Wl|| + ||W 2 || + (>Vl,W2) + ||e|| 

which can be substituted into 

? _ N-M b Q a -Q b 

M b -Ma Q b 

to yield 

^_ n — M b {\\w 2 \\ + {n -iy- 



(33) 



M b -M a 



Qi, 



- 1 



(34) 

where note that the expression depends 
solely on the signal generated by the planet 
with the weakest contribution. Its form is 
analogous to the one-planet case, where the 
modified variance a' 2 is defined as 

a 2 j_ (Wi,W 2 ) , WhW 
^ + 7a73TT + (N^T) ' 

. (35) 

In principle, the two additional terms in the 
right-hand side cannot be set to zero. First, 
the contribution from the two planets to the 
signal cannot be separated and when ||W || 
is obtained by means of an one-planet fit, 
the result is certainly contaminated by the 
fact that there is a second planet contribut- 
ing to the signal ||V||. Second, the covari- 
ance (WijW^) may be different from zero, 
especially if the two planets are assumed to 
be in resonance. However, both terms are 
inversely proportional to (N — 1), and it is 



expect that a' 2 
data sets. 



a 2 for sufficiently large 
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